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We develop a general numerical method to study the zero temperature properties of strongly 
correlated electron models on large lattices. The technique, which resembles Green's Function Monte 
Carlo, projects the ground state component from a trial wave function with no approximations. 
We use this method to determine the phase diagram of the two-dimensional t-J model, using the 
Maxwell construction to investigate electronic phase separation. The shell effects of fermions on 
finite-sized periodic lattices are minimized by keeping the number of electrons fixed at a closed-shell 
configuration and varying the size of the lattice. Results obtained for various electron numbers 
corresponding to different closed-shells indicate that the finite-size effects in our calculation are 
small. For any value of interaction strength, we find that there is always a value of the electron 
density above which the system can lower its energy by forming a two-component phase separated 
state. Our results are compared with other calculations on the t-J model. We find that the most 
accurate results are consistent with phase separation at all interaction strengths. 



PACS numbers: 71.10,+x, 71.45. Gm, 74.20.-z 
I. INTRODUCTION 

Correlated quantum many-body systems have pro- 
vided a host of new phenomena such as new states of 
matter, new forms of ordering transitions, etc. These 
phenomena are the result of the appearance of funda- 
mentally new relevant degrees of freedom which emerge 
as a coherent superposition of the underlying degrees of 
freedom of the many variable system. Once one identi- 
fies the important degrees of freedom, these degrees of 
freedom can be treated by analytical or semi-analytical 
techniques which are variants of generalized perturbation 
expansions around the defining framework. The problem 
which arises is that such frameworks cannot be imagined 
before hand unless there are hints from either experi- 
ment or numerical studies of models correctly capturing 
the dynamics of more basic degrees of freedom which, at 
first sight, seem featureless. 

Several such models on a discrete lattice exist and a va=. 
riety of numerical techniques are at our disposal to use.u 
Exact diagonalization techniques suffer from the fact that 
the dimensionality Nh of the Hilbert space grows expo- 
nentially with system size N s (number of sites). Taking 
into account all the symmetries of the problem can re- 
duce the size of the invariant subspaces to smaller size 
Nr (which may be a few orders of magnitude smaller 
than Nh)- However, the largest possible size increases 
only with the logarithm of the ratio of Nh/Nr. In par- 
ticular, most interesting quantities scale with the linear 
dimension of the system which scales with [Iti(Nh /Nr)} 3 
where d is the dimensionality of the problem. Renor- 
malization group approaches, such as the density matrix 
renormalization group (DMRG) technique^ have been 
very successful in one dimension but there are significant 
limitations in higher dimensions. 

Attractive alternatives seem to be stochastic methods 



such as quantum Monte Carlo which can give information 
about larger size systems. Many interesting problems, 
however, involve fermionic degrees of freedom. If one at- 
tempts a simulation of fermions at low temperature one 
encounters the so-called fermion sign problem. Namely, 
one needs to define configurations which carry a phase 
(a positive or negative sign) along with their statistical 
weights, a reflection of the transformation property of the 
fermion wave function under particle permutation. In the 
computation of many quantities of interest, such as the 
energy, the "positive" and the "negative" configurations 
give nearly opposite contributions, leading to wildly fluc- 
tuating weights. The negative-sign problem causes the 
statistical fluctuations to diverge exponentially with in- 
creasing system size for fixed density. 

The Green's function Monte Carlo (GFMC) method 
has been successfully applied to lattice spin systems, 
in particular to the square, lattice spin- 1/2 Heisenberg 
quantum antiferromagnet.HtjLl In this case, through 
the Marshall-sign transformation, the problem can be 
mapped to a hard-core boson problem which presents no 
sign problem and solved accurately on large size systems. 

An approximate method to deal with the sign prob- 
lem in fermionip systems is the fixed node (FN) 
approximationoa This approach projects a trial state 
onto the best variational state with the same nodal struc- 
ture, thus controlling the statistical fluctuations. The FN 
appi'ieximation has been used for lattice fermion systems 
alsoMB 

The GFMC method for lattice fermions without the 
fixed node approximation has been applied to one- 
dimension systems .E3 In two dimensions it has been ap- 
plied to the t-J modeLLn the limits of small numbers of 
electronaij or holes. E£E3 

In this paper, we present an efficient implementation 
of GFMC for lattice fermions at arbitrary densities of 
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electrons or holes. We demonstrate the utility of the 
method in the case of the two dimensional t- J model. The 
method projects a trial wave function onto the lowest en- 
ergy eigenstate that it overlaps. If the trial state overlaps 
the ground state, the projection yields the ground state. 
The projection becomes statistically more accurate as the 
ground-state component of the trial state increases rel- 
ative to the excited-state components. Results obtained 
for the t-J modeLwith this method have been published 
by the authors .Ell In this paper we present the general 
method and in addition, our results for the t-J model 
are presented in detail and compared with other recent 
calculations. 

The t-J model is thought to contain some impor- 
tant aspects of the environment in the copper-oxide 
superconductors. For instance, the calculated single- 
hole spectrumlifl is in— agreement with the results of the 
photo-emission data.Ea The model gives rise to a two- 
hole bound stateE£l with the d x i_ y i symmetry which is 
the believed symmetry of the superconducting state in 
these materials. In addition, Emery, Kivelson, and Lin 
(EKL)El suggested that the cuprates are near an elec- 
tronic phase separation (PS) instability which is pre- 
vented by the long-range part of the Coulomb interac- 
tion. In the phase-separated state, the holes cluster to- 
gether with a certain density of electrons, leaving the 
rest of the system in an antiferromagnetic state with no 
holes. Phase separation in the t-J model has been stud- 
ied by a number of techniques, which seem to be giv- 
ing conflicting .conclusions QEj l3 High temperature se» 



ries expansionscJi£j and some studies on small systemsc£l 
indicated that phase separation does not occur at the 
physical region for the cuprates, namely._J.Zi. ~ 0.3 — 0.4, 
while other studies on small systemso'EII found this 
region was unstable to phase separation. Using the 
GFMC method presented in this paper, Hellberg and 
ManousakisEH found that the t-J model has a region of 
phase separation at all interaction strengths, .—. 

In recent work, Calandra, Becca, and SorellaEH (CBS) 
emphasize that the phase separated region does not ex- 
tend below J/t < 0.4. In addition, the DMRG method 
has been also acplied to this problem by White and 
Scalapino (WS)S__I who find that the ground state of 
the t-J model on the square lattice is characterized by 
stripes. 

In the last section of this paper, we compare our re- 
sults to those of CBS and WS, concluding that the phys- 
ical region of the model is very close to or inside of a PS 
instability. Namely, the early conclusions that the phys- 
ical region of the t-J model is far from the critical J c for 
phase separation, are largely invalid. 



II. NUMERICAL METHOD 

Even though our formalism is general and can be ap- 
plied to other lattice fermion models, we shall use the 



example of the t-J Hamiltonian on a two-dimensional 
square lattice. This model in itself is a non-trivial ex-, 
tension of the square lattice Heisenberg antiferromagnetlil 
where GFMC was first applied on a lattice model. 

The t-J model is written in the subspace with no dou- 
bly occupied sites as 
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Here (ij) enumerates neighboring sites on a square lat- 
tice, c\ a creates an electron of spin a on site i, n^, = 
E C r c ler c io-) an d Si is the spin-i operator. 



A. Details of the Projection 

We take a trial wave function 'J and project it onto the 
ground state by generating a series of increasingly accu- 
rate approximants to the ground state labeled by integers 
|m) = (H - W) m |$). Here H is the Hamiltonian and W 
is an appropriately chosen numerical constant .t__~__l 

We may expand the trial state in terms of the exact 
eigenstates: 
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where |$o) is the ground state, |*i) is the first excited 
state, and the Oj's are the expansion coefficients. Rewrit- 
ing the projected states in this way, we see 

\m) = (H - W) m \V) (3) 
= a Q (E Q - W) m \$ ) + ax{E x - W) m \$i) + ■ ■ • (4) 
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where Ei is the energy of the i'th eigenstate. So |m) 
approaches the ground state for large m provided 
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for all excited state energies £?i>o- The projection can 
be formulated in this simple way because eigenvalues of 
lattice Hamiltonians are bounded from below and above. 
Continuum problems require a different form for the pro- 
jection operator Ptt 5 ! In what follows, we assume the off- 
set constant W is incorporated in the Hamiltonian. 

From (|^) we see rate of convergence with m is governed 
by the overlap of the trial state with the ground state 
and the energy of the lowest excited state overlapping 
the trial state. In Section II D, we describe the steps 



taken to insure fast convergence. 

To calculate ground state expectation values of an ar- 
bitrary operator A, we take the large m limit of 
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For large values of m, we cannot evaluate H m directly. 
The number of position-space states generated diverges 
exponentially with the power m, so we calculate H m by 
a stochastic method similar to Neumann-Ulam matrix 
inversion. c3 We decompose H into a product of a transi- 
tion probability p a p to make a transition from state a to 
state P and a residual weight w a p as 



where 
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stochastically, we average over m-step random walks 
«o — > o;i — > • • ■ — > a m _i — > a m , where each a, is a 
position space state, giving each walk the accumulated 
weight 

(11) 

The probability of the walk ao — > ai, • ■ ■ — ► a m is 

P(.OtO ; ) <^m ) — Pao&iPaia2 ' ' ' Pcx m — ia m • (1^) 

Thus, it follows that 

(aol-H^lom) = J! W{a ,ax,...,a m ) (13) 

for a large number of walks guided by the probability 
B. Importance Sampling 



The Monte Carlo sum (|13|) is evaluated most efficiently 
using importance sampling. We cannot use the trial state 
as a guiding function for the random walk, since the guid- 
ing function must be positive for all allowed states. La- 
beling our guiding function \I> , we let 
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where the normalization is simply the local energy: 
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Defined in this way, (|TJ) satisfies (g), and the residual 
weight is 



(16) 



resulting in the accumulated weight for the m-step walk 
given by Eq. Q. 

For an antisymmetric trial wave function 4 ,T , the stan- 
dard algorithm to evaluate 
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where vp^* = (^> T \a), is to generate a set of M ini- 
tial states {ai} with probabilities proportional to Q ai oc 
| 2 using Metropolis sampling as in Variational Monte 
Carlo. At each initial state \oti), we start an m-step ran- 
dom walk, ending in the state For large M, 
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C. A New Approach 

The standard algorithm is inefficient since a random 
walk in configuration space of length m must be gener- 
ated for each term in the sum ([l^). The details of the 
intermediate states are thrown away. The expectation 
value (^\H m \^) can be evaluated more efficiently if the 
generation of the initial states {ai} is combined with the 
generation of the random walks. In the random walk, 
new states are chosen with a probability given by (|H|). 
After a large number of steps, these states are distributed 
with a probability 
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which is derived by solving the "detailed balance" condi- 
tion 
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where p a p is the probability to make a transition from 
configuration a to /3 given h« (|l4|), and Q a is the prob- 
ability to visit a state |a)£3 Thus we may use states 
generated in the m-step random walk as initial states for 
new m-step random walks. 

For maximum efficiency, we use every state generated 
as the starting point for a new walk, so at each step we 
calculate different stages of m-walks simultaneously. We 
simply generate one very long random walk using the 
probability (|l4|). At each step in the walk, we look m 
steps into the past to evaluate an element of (|l^). The 
computer time needed to calculate a given number of 
observations of (H m ) is independent of m. An additional 
advantage is that since only one long random walk is 
generated, we may calculate all different powers m in 
parallel. The fundamental observation becomes 
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The method is easily generalized to evaluate the ex- 
pectation value A m = (^\H m AH m \^) for any diagonal 
operator A, such as the density or spin structure factors, 
riiUj and SfSj. At each stage in the walk, we look m 
steps into the past to obtain the expectation value of (A) 
and 2m steps into the past to calculate the accumulated 
weight. By summing M observations from a walk, we 
find 
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where j = i — m and k = i — 2m. 

The speed of convergence of the procedure with power 
m is determined by the ratio R — \Ex — W\/\E — W\ < 1, 
where E\ is the energy of the first excited state overlap- 
ping the trial state. Since this gap is caused by the finite 
size of the system, we generally calculated powers of the 
Hamiltonian up to several times the linear system size. 



D. Trial Wave Functions 

Care is needed to choose a trial state with maximal 
overlap with the true ground state. We restrict ourselves 
to total spin singlet states with zero momentum and try 
to write a very arbitrary form yielding a good initial guess 
throughout the phase diagram. 

We use a Jastrow resonating- valence-bond wave func- 
tion for the trial state, written 

II /(r^-r^OlRVB) (23) 

i<J,0",fT y 

i<j,er,cr' k 

where c k(T is the usual Fermion creation operator and Pn 
projects the state onto the subspace with the number of 
particles fixed to be N. 

It is important that the Jastrow factor / correlate all 
pairs of particles independent of spin, yielding a corre- 
lated state that is still a total spin singlet. If we allow dif- 
ferent Jastrow factors for like and unlike spins, we could 
usually reduce the variational energy of the trial state. 
However, the resultant state would be a superposition of 
many spin states and in general would overlap excited 
states with non-zero spin closer in energy to the ground 
state than the lowest spin zero excited state, resulting in 
slower projection than a singlet trial state with higher 
variational energy. 

We write the_determinantal part of the trial state in 
the usual way £3 The ratio = t'k/«k is the physical 
quantity, and, assuming ak = a_k, we define a(r) as its 
Fourier transform: 



<z(r) = at cos(k ■ r). 

k 

Then 

|RVB) = P N []( Uk + v k c^ci kl )\0) 
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can be written as the — x y determinant 
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in position space. 

In this form, |RVB) spans a broad class of Fermion 
wave functions. A Fermi liquid state corresponds to 



Iks Fermi sea 
otherwise 



(27) 



while by allowing other choices for ak, the wave function 
can describe a pairing state, which may be s-wave, d- 
wave, or something more general. 



E. Guiding Function 

We are tempted to use the magnitude of the trial state 
as our guiding wave function, but this would be a serious 
mistake. By construction, the sites of a periodic lattice 
lie at high symmetry points, and the nodes of a Fer- 
mion wave function also respect these symmetries. One 
finds (Nf'Jla) = for a significant fraction of states in the 
Hilbert space not violating the Pauli exclusion principle. 
Since 

(* t \hv\* t ) = J2 (C |ffK>-->™-i|H|Cj 

(28) 

where the intermediate sums over ai, a% . . . , a n -i span 
the complete Hilbert space, we must guide with a positive 
function. 

Every guiding function that samples the complete 
Hilbert space will yield correct results. Our challenge 
is to pick a function that minimizes the statistical fluctu- 
ations of our output. The guiding function cannot ame- 
liorate the sign problem in (|2l|). We can, however, choose 
a function to reduce the fluctuations in |$ T |/|W G |. We 
define 
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FIG. 1. Schematic behavior of the guiding function near 
a node. The squares and diamonds are the trial state and 
its negative, respectively. The circles are the bosonic state, 
and the guiding function, VP = max {|* T |,cvf s }, is shown 
by the filled symbols, r represents a single coordinate of one 
electron. 



* G = max{|* T |,c* B } 



(29) 



where VP is a positive function, typically a good varia- 
tional state of the bosonic Hamiltonian— We take \E' B to 
be a spin-dependent Jastrow function.E-3 This is similar 
to a choice used in continuum problems, but on a discrete 
lattice, it is not necessary to match the first derivatives 
as in the continuum.La 

We rescale c so the effective number of configurations 
contributing to the norm is approximately N ps l/L 
for an L x L system.EJ This guiding function is shown 
schematically in Fig. |l| 

For the guiding function, we use the Jastrow-pairing 
function 
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of Bose spin i particles (i.e., two kinds of bosons, up 
bosons and down bosons) . We have chosen this function 
to mimic the physics of the Fermion state as closely as 
possible without having nodes. Since it is not impor- 
tant to guide with a spin singlet function, we use a more 
arbitrary spin-dependent Jastrow factor where like spin 
particles are correlated differently than opposite spin par- 
ticles. 



F. Walking Through The Nodes 



the determinant and inverse may be updated after such 
moves efficiently in 0(N 2 ) steps. This so called "inverse 
update" works well for variation Monte Carlo or fixed 
node Monte Carlo, but it cannot be used directly with 
GFMC on a lattice since the random walk steps directly 
on nodes for a significant fraction of steps. In a reason- 
ably dense system, we find as many as 1/3 of the steps 
land on nodes. On a node, the matrix is singular, the de- 
terminant is zero, and the inverse is undefined. Recalcu- 
lating the determinant and inverse after walking through 
a node will cause the running time to scale as 0(N 3 ). 

We developed a new 0(N 2 ) technique to hop over 
nodes without recalculation of the determinant or in- 
verse. The essence of the method is this: When the ran- 
dom walk generated by the guiding function hits a state 
or series of states where the determinant of the trial func- 
tion vanishes, we generate a "detour" walk around the 
region where the matrix is singular, rejoining the guid- 
ing walk when the determinant is non-zero again. We 
stress that the real random walk goes though the node, 
the detour walk is a fictitious one which is used only to 
calculate the determinant and its inverse. It serves only 
as a calculational tool for the inverse update. The details 
of this detour-walk-approach of evaluating the determi- 
nant and its inverse when the walk went through a node 
are explained in the Appendix |a|. 

Since we often land on a node where the inverse of the 
matrix ( ^6| ) is not defined, it is difficult to calculate the 
probabilities p a p (Eqs. [l4| and |l5|) for the random walk 
we defined earlier. When we are on a node with iff^ = 0, 
we need to choose the next step of the walk from the 
various possible (3 states. Therefore we need to calculate 
z a with Eq. (|l5|). Since we are on the node, calculating 
each requires a detour walk which in itself takes TV 2 
steps. Thus, the calculation of z a is an 0(N 3 ) process, 
whereas when ty^ ^ it is an 0(N 2 ) process. Therefore, 
we make the following adjustments so that each step is 
an 0(N 2 ) process. We define 
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(31) 
(32) 



(33) 



and if = 0, 



To evaluate the determinantal wave function, at the 
beginning of the random walk we calculate the deter- 
minant, an 0(N 3 ) operation for a N x N determinant, 
and its inverse, also an 0(N 3 ) operation. Each kinetic 
step in the walk changes either a row or column of the 
determinant, while superexchange changes, both a row 
and column. Ceperley, Chester, and Kalosci showed that 



fa/3 



(34) 



It is easy to show that detailed balance is obeyed by 
these definitions. The advantage of using these probabil- 
ities is that if = 0, then calculating 'J/T for all (3 with 
H a f3 ^ not required. Since we don't have the inverse 
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matrix of the trial function's determinant in state a, such 
a calculation would be computationally very expensive. 

Calculating , 3> ( f is always easy due to the simple form 
of (|30|), but is more difficult. The kinetic operator 
moves a single electron, so Vf^ may be calculated in O(N) 
steps since the inverse need not be updated. The superex- 
change operator moves two electrons and changes both a 
row and column of (^6j). Updating the inverse to calcu- 
late this term requires 0{N 2 ) steps, which would cause 
the overall algorithm to scale as 0(N 3 ) for the system. 
In Appendix [b| we derive an O(N) method of calculating 
superexchange. 



G. Trial State Optimization 

It is important that we start the GFMC with good 
trial and guiding states. In this section, we describe our 
method for optimizing these functions. 

In continuum systems, one usually assumes a func- 
tional form for the trial and guiding functions and opti- 
mizes a function of the energy to find the best variational 
parameters. On a lattice, there are only a finite number 
of distances r or equivalently wave vectors k in any given 
simulation, so we allow the functions in ( ^3| ) and (PQ ) to 
have a parameter describing each distance or wave vector 
not related by symmetry. 

For the Jastrow and position space pair factors, /(r) 
and g(r), we apply all rotational and mirror symmetries. 
Translational symmetry is always assumed. However, we 
insist only on the mirror symmetries about the axes for 
the Fermion pairing field at, so the function may be any 
linear combination of an s and d x 2_ y 2 pairing state. The 
mirror symmetry excludes d xy symmetry. 

For a 20 x 20 lattice, we have 172 parameters for the 
trial state and 192 parameters for the guiding function 
with the pairing term. 

We tried optimizing several functions of the energy, 
but found minimizingj-the variance of the local energy 
to be the most robust .c3 We generate a set of configura- 
tions {ai, «2, ■ ■ ■ ! ct m } distributed according to a weight 
w ai . The configurations remain the same throughout the 
minimization procedure. We minimize the function 
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where E is a guess for the ground state energy that we 
determine self consistently. We use the same function to 
optimize both our trial and guiding functions. 

With a finite random walk, the calculation of the en- 
ergy in (|3^) uses many more states than the calculation of 
the norm. Occasionally, this created instabilities, which 
we cured by deriving another way of calculating the norm 
using all the neighbors in the random walk. We may write 



= £ [\* a \ 2 (i-A a ) + B a J2 l^l 2 ) > :!7i 

where {H a } is the set of all states neighboring \a) by 
application of the Hamiltonian. We see ( p7|) follows from 
( |36| ) if we choose B a = C and A a = CN a for some 
constant C, where N a is the number of neighbors of \a) 
where 'J does not vanish. Since this version of the norm 
is calculated from all the states entering the energy, no 
factors in the numerator of (35) are absent from the de- 
nominator. 

We calculate the effective number of configurations 
contributing to the normalization as 




(38) 



This quantity approaches n if all states contribute equally 
to the integral and drops to 1 as one state begins to 
dominate.cJ We adjust the length of our random walks 
so JV e ff is at least 10 times the number of parameters 
being optimized. 

We found that close to phase separation, the standard 
Metropolis algorithm develops a small acceptance ratio, 
and tends to stay in the same configuration for many 
steps. In order to sample more phase space quickly, we 
choose our configurations using the transition probabil- 
ity (|l4|) where we take H to be the off-diagonal part of 
the Hamiltonian, ensuring a new configuration with each 
move. Thus the configurations are distributed according 
to the weight w ai = z a . |*g 4 1 2 , where z ai is given by @. 



III. FITTING PROJECTION OUTPUT: INVERSE 
THEORY 

A. The Ground State Energy 

The Green's Function Monte Carlo procedure takes a 
trial state and projects it onto the exact ground state. 
Its output consists of the observables for the energy 



#(") = mH n m 
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where the trial state |\&) has been normalized. For any 
operator A which does not commute with the Hamilto- 
nian, using the present Monte Carlo method we calculate 



mH m AH n m 
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as functions of powers of the Hamiltonian m and n. 

These values converge to their ground state values, ex- 
cept for a normalization factor, for large powers n and m. 
However, their statistical errors increase exponentially 
with increasing power due to the Fermion sign problem. 

To extract the most information on the ground state, 
we use the calculated observable for all powers less than 
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some maximum power p ma x- By including the highly 
converged small powers in the approach, we obtain much 
more accurate ground state properties than can be ob- 
tained using the large powers alone. Let us consider the 
ground state energy as an example to demonstrate the 
approach next. 

Let us define the T-spectral function c(£) with respect 
to the trial state \^) as 
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To show this, one may expand the trial state in the exact 
eigenstates of H as 



(42) 



It is immediately evident from the above that the poles 
of c(£ ) and of the exact spectral function are at the same 
energy values for those eigenstates which have non- 
zero overlap with the trial state |\&). 

In order to proceed we discretize the energy interval us- 
ing a fine mesh with AE, thus, the energy takes discrete 
values E* m , m = 1, 2, M and the T-spectral function 
c*(E) is written as 
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where c* m > are non-negative real numbers. Thus, this 
spectral function is thought of as a histogram, where in 
each fine slice of the histogram the value of the inte- 
gral of c(E) multiplied by any function f(E) is simply 
c* m f(E* m ). We have used a mesh interval AE smaller 
than the finite-size gap between the lowest and the first 
excited state. Thus, the contribution of the ground state 
to c* is accurately represented as a single delta-function 
peak. Namely, up to some m = mo, c* n<mo = 0, while 

C mo > and C ma<m<m 1 = and C m± > 0, etc. 

Then the moments of the T-spectral function c*(E) 
can be calculated using both Eq. (|4l]) and Eq. (Q) and, 
thus, we obtain: 



A I 



(9\H»\Hf) = E C ™ E * 



(44) 



m— 1 



where n = 0, 1, ...,p ma x- Since p m ax < M, because typi- 
cally p m ax = 40— 60 and M — 200 — 500, we have more 
unknowns (the c^'s) than equations. However, the solu- 
tion needs to satisfy the constraint c*„ > which limits 
the possible solutions. The optimal way to find the most 
likely solution is to minimize the x 2 . 

We gain very large computational savings by calculat- 
ing all powers of H in parallel. However, this results in 
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FIG. 2. The T-spectral function of the full Hamiltonian. 
The lowest energy peak is at the lowest energy eigenstate of 
the system which is non-orthogonal to the trial state. The 
value of spectral weight is the square of the overlap of the 
true ground state to the initial trial state. 



statistical correlations between results of different powers 
which must be treated accordingly.^ 

We divide the measurements into M bins. The covari- 
ance matrix is defined 



1 



1 



M - 1 



-E&) ) (45) 



fc=i 



where E^ is the average of the i'th power in the fc'th 
bin. For uncorrelated output, C is diagonal. With corre- 
lations, x 2 i s defined 



X 



= J2( E(l) - E (l) *)C^(E^ - E^*), 



(46) 



where E^* is the fitting function given in terms of the 
coefficients c m (which are to be determined by this min- 
imization) by means of Eq. (|44|). 

When C is diagonal, its inverse is trivial. For more 
general C, small errors in its components can result in 
large errors in its inverse, so it is important to calculate 
C accurately. Increasing the number of bins decreases the 
statistical error in C but increases the systematic error 
due to autocorrelations. To balance these two sources 
of error, we choose the number of measurements in each 
bin to be n = Afp m _. where p max is the maximum power 
of the Hamiltoniano We calculate statistical errors with 
the bootstrap methodJU 

Fig. H gives a typical example for the T-spectral func- 
tion c*(E) obtained from the calculation of the two- 
dimensional t-J model. The lowest value of E* where 
we have a delta function peak gives the lowest eigenstate 
of H which is not orthogonal to the trial state. The value 
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FIG. 3. The energy for 10 electrons in a 4 x 4 lattice as 
a function of the power of the Hamiltonian. The value of 
the energy obtained by the extrapolation method described 
in this section is also shown. 



of the peak gives the square of the overlap of the lowest 
energy state to the trial state. 

We have tested our method by comparing our results 
with exact results for the 4x4 size lattice with several 
electrons. In Fig. || we show the results for the energy 
as a function of the iteration for the case of 10 electrons. 

The energy estimate defined by 



En 



(*| J ffP+ 1 |*) 



(47) 



as a function of p is shown in Fig. (|3|) starting from 
the Jastrow-RVB value of the energy at m = 0. Notice 
that with the length of our walk in configuration space 
which we used for this calculation, the error (which al- 
ways grows exponentially) becomes annoyingly big for 
values of p not shown. The value of the energy obtained 
by the extrapolation method described in previous sec- 
tion is also shown. By using the information contained 
in all the powers of H up to p m ax instead of the just the 
estimates of the energy at or just below p ma x we obtain 
a much better estimate for the energy. 



B. Other Operators 

For an arbitrary operator A, we have 
(*|iTAff«|*) = Y.iEiE^m^mA^i)^^) 

ij 

= J dEE 2p a(E) (48) 

where the operator overlap function a(E) is given by 
a(E)=J2S(E+^E7E~)(^ i )(^ i \A\^ j )(^). (49) 
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FIG. 4. The energy for 50 electrons in a 7 x 8 lattice as 
a function of the power of the Hamiltonian. The value of 
the energy obtained by the extrapolation method described 
in previous section is also shown. 



Following the approach for the energy T-spectral func- 
tion, we can define a discrete overlap function 



M 



a*(E)=J2S(E-E*)a*. 



(50) 



i=0 



Here the values of E where the T-spectral function a*(E) 
attains peaks are all possible geometric means U EiEj 
of all the eigenenergies which correspond to eigenstates 
which have non-zero overlap with | , F) and they give non- 
zero matrix element of A. The lowest energy peak corre- 
sponds to the geometric mean of the ground state energy 
with itself, i.e Eq = ^/EqEq thus, if the ground state is 
not degenerate it is uniquely specified. Here we also need 
to solve for all the a* given that they obey the following 
Pmax equations: 

A I 



i=0 
M 

(*\HAH\V) = 



i=0 



A I 



(V\H p AH p \iB) = J2(E*) 



2/V 



(51) 



Fig. g gives a typical example for the spectral 
function a*(E) obtained for the spin-structure function 
S(tt/2,tt/2) from the calculation of the two-dimensional 
the t-J model. Notice that the energy of the lowest peak 
is the energy of the lowest energy state having non-zero 
overlap with the trial state and which has non-zero ma- 
trix elements with the operator A. For the value of the 
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FIG. 5. The T-spectral function associated with the 
spin-structure function. Notice that the lowest energy peak is 
at the same energy as that of the energy T-spectral function. 
The value of spectral weight is related to the spin-structure 
function in a simple way. 

peak, which is a* , the expectation value of A can be cal- 
culated in a straightforward manner. 

In principle, wc can also extract information about the 
excited states along with the ground state. This possibil- 
ity is indicated by the fact that we can see higher energy 
peaks in these spectral functions. 

In Fig. § we compare our results for the spin-spin cor- 
relation function with that obtained by exact diagonal- 
ization of the 4x4 size system with 10 electrons. Notice 
the fine scale used to be able to distinguish the difference 
between exact and extrapolated correlation functions. 

IV. RESULTS FOR THE 2D T-J MODEL 

A. Maxwell Construction: Phase Separation 

A number of methods to determine the phase sepa- 
ration boundary numerically have been used. In one di- 
mensional systems the divergence of the density structure 
factor at long wavelengths has been used.E3 Divergence of 
the compressibility as determined from the second deriva- 
tive of the energy with respect to electron or hole density 
was used successfully on the one-dimensional t-J model 
by calculating the energy for three differing densities. c£l 
In the one-dimensional model, phase separation occurs 
between two regimes, one with no electrons while the 
other contains some electrons and some holes. For a fi- 
nite system, electrons may tunnel through the vacuum, 
lowering the ground state energy. For this reason, the 
inverse compressibility actually passes through zero and 
becomes slightly negative. This effect is a surface effect 
and vanishes in the limit of infinite system size. 
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FIG. 6. The spin-spin correlation function for 10 electrons 
in a 4 x 4 lattice as a function of the power of the Hamiltonian. 
The value of the correlation obtained by the extrapolation 
method described in previous section is also shown. 

In the one-dimensional t-J model, the compressibility 
diverges continuously at the transition point in contrast, 
to the discontinuous transition in two dimensions O'H 
The behavior of the energy derivatives across the phase 
separation boundary is discussed in Appendix |c[ We 
have verified that in one dimension, the Maxwell con- 
struction yields the same phase-separation boundary as 
that calculated using the inverse compressibility. 

In the two-dimensional t-J model, the situation is more 
complicated. The Fermi surface can change dramati- 
cally with electron density for a given system size. These 
strong shell effects make accurate comparisons of energies 
calculated with different numbers of electrons impossible. 

Many of the previous studies used a vanishing inverse 
compressibility as the criterion for the onset of phase 
separation .oE3 The compressibility, however, is not the 
proper observable to find the phase-separation boundary 
in the two-dimensional t-J model, where the transition is 
first order. It is true that the compressibility diverges in 
the region of phase separation, but it jumps discontinu- 
ously at the boundary with the uniform phase. Numer- 
ically, this discontinuity is difficult to see in even large 
finite systems due to the surface energy of the two co- 
existing phases. When one is in the region of the phase 
diagram where phase separation exists, then, the com- 
pressibility suffers strong finite-size effects because of the 
rather large surface energy of the two coexisting phases. 

The ground-state energy as a function of electron den- 
sity at J = 2.5t for 32 electrons on a variety of system 
sizes is shown in Fig. |7]. These finite systems necessarily 
constrain the electron density to be uniform on the length 
scales of the system size. We fit the discrete data to a 
polynomial, e(n e ), shown as the solid curve, in order to 
treat the energy as a continuous function of density. The 
dashed line, e ps (n e ), is a linear function that intersects 
the Heisenberg energy, en at electron density n e = 1 and 
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intersects e(n e ) tangentially at a density labeled n ps . 

It is straight forward to show that the ground state of 
the infinite system at a density n e > n ps cannot be a 
uniform phase, because the energy of the uniform phase, 
e(n e ), is higher than e ps (n e ) at the same density. This 
latter energy corresponds to the energy of a mixture of 
two phases, one at electron density ua = 1 and the other 
at electron density ns — n ps . Therefore the infinite sys- 
tem phase separates into two regions with densities ua 
and riB, and its ground-state energy is given by e ps {n e ), 
the value of the dashed line at the average density of-4hc 
system. This is known as the Maxwell construction.E3 
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FIG. 7. The ground-state energy per site at J = 2.5t 
for 32 electrons. For clarity, the energies are shifted by a 
linear factor, — enn e . The circles with error bars show the 
energies calculated on lattices of dimensions 6 x 6, 7 x 6, 
28 x 28. A sixth-order polynomial fit to the data is shown 
as the solid line, which is extended to the Heisenberg energy, 
the square at energy zero in this shifted plot. The dashed 
line shows the ground-state energy of the infinite system in the 
phase-separated region. We find the onset of phase separation 
occurs at n ps = 0.296±0.004, while the inverse compressibility 



vanishes at n co „ 



0.52 ±0.10. 



The energy of the infinite system is given by the solid 
line in Fig. |7 for n e < n ps and by the dashed line for 
n e > n ps . A difference between the Maxwell construction 
in Fig. [7| and that commonly used is that the density of 
one of the constituent phases, the Heisenberg phase at 
n e = 1) lies at an extreme limit of the allowed density 
range. It is not possible to add electrons to the Heisen- 
berg solid, which has one electron on every site, so the 
dashed line is not tangent to the fitting curve at n e = n psi 
it is not at n e = 1. If the t- J model did allow electron den- 
sities n e > 1, then the intersection point of the solid and 
dashed lines would be shifted to higher densities where 
the curves could intersect tangentially. The dashed line 
might intersect the solid curve tangentially at the higher 
density point in this region. At any electron density in 



ergy from that of the uniform phase approximated by 
the fitting polynomial in Fig. 0, by separating into two 



regions with densities ua 



< ps 



and tib = 1, resulting in 



the range, 



<'pS 



< n e < 1, the system can reduce its en- 



an energy given by the dashed line at the average density. 

In order to be stable, the energy of the infinite system 
must be concave everywhere. Given the solid line in Fig. 
fj] and the allowed density range of the t-J model, the 
dashed line drawn in the figure is the only line possible to 
make the energy of the infinite system globally concave. 
This energy is given by the solid line for n e < n ps and 
the dashed line for n e > n ps . 

Since the t- Jmodel is denned only for electron densities 
> n e > 1, the energy of the phase-separated state, the 
dashed line in Fig. [?], is tangent to the fitting curve at 
n e = n ps but not at the Heisenberg point, n e = 1, an 
extremal density. 

We never examined systems with densities n e > 0.94, 
so we cannot exclude the reentrance of a homogeneous 
phase in this region. For such a phase to be stabilized, 
the solid curve in Fig. [j] would have to drop back below 
the dashed line in this density range. We never saw any 
indication of this possibility at any interaction strength. 
If a reentrant homogeneous phase did exist, the phase- 
separated region would persist at densities n e < 0.94. 
The new Maxwell line would lie slightly below the one 
drawn and would be tangent to the solid curve at both 
intersections, but the phase-separated region would per- 
sist at densities n e < 0.94. We never saw any indication 
of this possibility at any interaction strength. 

In the infinite system, the energy is a concave func- 
tion of density with continuous first derivative. In the 
phase-separated regime, the energy linearly interpolates 
between the energies of the two constituent phases. The 
inverse compressibility, the second derivative of the en- 
ergy with respect to density, is positive in the uniform 
phase and zero where the system phase separates, in 
the phase-separated regime. On a finite lattice, the sur- 
face effects of phase separation raise the energy, and we 
find the inverse compressibility calculated from differ- 
ent system sizes obtained from a fit to the discrete data 
remains positive even where the system is unstable to 
phase separation. For electron density n e , the completely 
phase-separated energy per site, e var = enrie, is a vari- 
ational upper bound to the phase-separated energy, is 
given by the completely phase-separated energy per site, 
£var — znn e . To show that, let N e be the total number of 
electrons on a square lattice of N sites (n e = N e /N < 1). 
Then, the total energy of the system is less or equal to 
e/r-^e plus corrections which are negligible compared to 
this term in the thermodynamic limit. To show that, let 
one constrain all the electrons to be close-packed in one 
region of the lattice so that only the interaction term in 
(0) is operative. If the system chooses not to do that, it 
means that the energy of any other chosen state is lower 
than the energy of this artificially imposed state. Thus, 
the energy per site of the unconstrained system is lower 
than eu\n e . The system is definitely unstable to phase 
separation at least where this variational energy is lower 
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than the calculated energy per site of a finite system. 

We fit the calculated energies in the uniform phase to a 
low-order polynomial. The energies in the uniform phase, 
where there is no surface energy, suffer much smaller 
finite-size effects than energies in the phase-separated 
regime. We extrapolate the tangent of the fitting func- 
tion to electron density n e = 1. At the phase-separation 
boundary, the tangent to the curve e(n e ) extrapolates to 
the Heisenberg energy, e(n e = 1) = e#. This construc- 
tion ensures that the resulting infinite-system energy is 
concave everywhere. This is the only possible construc- 
tion that ensures the resulting infinite-system energy is 
concave everywhere. An example for J = 2.5< is shown 
in Fig. @. 



B. Results at J = t 
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FIG. 8. The ground-state energy per site at J = t for 
N e = 42 electrons and N s = 49, 56, 64, 72, 81, 90 sites (bottom 
curve) N e = 50 and N s = 56,64,72,81,90 (second from the 
bottom), and N e = 60 and N s = 64,72,81,90,100,110,121 
sites (top curve). Each curve has been shifted upwards with 
respect to the previous by 0.025 in order to distinguish them. 

Using the Maxwell construction we have determined 
the boundary for phase separation of the t- J model in the 
J < t region. In these calculations we minimize the shell- 
effects when varying the electron-density by keeping the 
number of electrons fixed at a closed shell configuration 
and changing the size of the lattice. In addition, we also 
choose the number of electrons to correspond to a closed 
shell configuration. This choice eliminates possible de- 
generacies of states at the Fermi level. Such degeneracies 
might favor flatness of the energy as a function of density 
which might be mistaken for phase separation. When one 
varies the number of electrons keeping the system size 
fixed when adding an electron to a closed shell configura- 
tion the kinetic energy goes up by a finite amount which 
leads to oscillatory behavior of the energy per site versus 
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FIG. 9. The energy per hole at J = t for N e = 42 (open 
diamonds) N e — 50 (solid squares) and N e — 60 (solid circles). 
The shell-effects have non-monotonic influence on the scaling 
with size of the energy per hole. 

density. An additional (technical reason) for wanting to 
keep closed shell configurations is that on a finite lattice 
it leads to an energy gap between the ground state and 
the first excited state which helps our projection method 
to converge. 

In Fig. U we study the size dependence of the critical 
value of electron density n ps for phase separation. The 
value of n ps is determined by making a cubic polynomial 
fit to each curve and using the corresponding energy for 
the lattice full of electrons calculated for the same num- 
ber of electrons and the same boundary conditions. This 
can be done by using the GFMC results obtained for vari- 
ous size lattices for the square- lattice spin-l/2-,Heisenberg 
antiferromagnet and using the extrapolation^! 



E/N = e + AA^ 3/2 



(52) 



where N is the total number of electrons. 

The ground-state energy per site at J = t for fixed 
number of electrons N e and for various number of sites 
N s can be grouped together. The bottom curve in Fig. 
H gives the energy per site as a function of density n e = 
N e /N s for N e = 42 and N s = 49(7 x 7), 56(7 x 8), 64(8 x 
8), 72(8 x 9), 81(9 x 9) and 90(9 x 10). The second from 
the top gives the energy per site shifted by a constant 
amount of 0.025 for N e = 50 and N s = 56, 64, 72, 81, 90. 
The third from the top gives the energy per site shifted by 
an 0.05 for N e = 60 and N s = 64, 72, 81, 90, 100, 110, 121. 
The unshifted energies are given in Table |. If we plot 
these curves on the same scale without a shift, they all 
fall on almost the same curve. However, there are small 
deviations from such a curve which are shell effects. 

Also shown in Table || is the lowest energy obtained by 
diagonalizing within the subspace spanned by a subset of 
2, 3, or 4 different powers. This is similar to the Lanczos 
algorithm for Quantum Monte Carlo data introduced by 
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TABLE I. The energy per site for J/t — 1 and for various 
electrons densities and size lattices. The extrapolation en- 
ergy is obtained by the procedure described in section III-A. 
The Lanczos energy is the lowest obtained by diagonalization 
within a subspace of projection powers. 





N s 


n e 


EExtrap/N s 


ELanc/N s 


32 


100 


0.32 


-0.9019(12) 


-0.9008(7) 




81 


0.395 


-1.0204(34) 


-1.0197(47) 




72 


0.444 


-1.0767(12) 


-1.0750(15) 




64 


0.5 


-1.1249(48) 


-1.1220(27) 




56 


0.571 


-1.1727(48) 


-1.1760(72) 




49 


0.653 


-1.1828(44) 


-1.1840(35) 




42 


0.762 


-1.1853(57) 


-1.1801(18) 




36 


0.889 


-1.1637(40) 


-1.1681(29) 


42 


121 


0.347 


-0.9463(13) 


-0.9448(5) 




100 


0.42 


-1.0583(38) 


-1.0548(39) 




81 


0.519 


-1.1335(37) 


-1.1383(41) 




72 


0.583 


-1.1562(34) 


-1.1548(15) 




64 


0.656 


-1.1695(18) 


-1.1662(22) 




56 


0.75 


-1.1768(32) 


-1.1765(24) 




49 


0.857 


-1.1674(33) 


-1.1656(16) 


50 


121 


0.413 


-1.0229(36) 


-1.0192(8) 




100 


0.5 


-1.1080(38) 


-1.1022(35) 




90 


0.556 


-1.1408(49) 


-1.1338(32) 




81 


0.617 


-1.1566(39) 


-1.1612(63) 




72 


0.694 


-1.1766(58) 


-1.1741(37) 




64 


0.781 


-1.1787(89) 


-1.1663(38) 




56 


0.893 


-1.1633(38) 


-1.1606(30) 


60 


110 


0.545 


-1.1451(40) 


-1.1442(35) 




100 


0.6 


-1.1601(36) 


-1.1602(18) 




90 


0.667 


-1.1698(34) 


-1.1697(23) 




81 


0.741 


-1.1769(28) 


-1.1759(41) 




72 


0.833 


-1.1771(37) 


-1.1710(35) 




64 


0.938 


-1.1628(18) 


-1.1612(15) 



Caffarel, Gadea, and Ceperley.a The energies obtained 
are upper bounds to the ground-state energy. 

We can greatly eliminate the shell effects by examining 
several size lattices but keeping the number of electrons 
fixed. We first fit each curve generated for fixed N e with 
a quartic spline. Next from the point e(n e — 1) on the 
graph, we construct the tangent to the spline which fits 
our points. The value of the density at which the tangent 
and the spline meet gives us an estimate of the phase sep- 
aration density for the given value of J/t. Individually, 
each number of electrons is consistent with a solution 
where a line beginning from the no- hole limit (n e = 1) 
and being tangent on the polynomial (that fits the data 
points) near n e — 0.745 (see Fig. ||). Therefore we con- 
clude that the finite-size effects are small in our method 
of determining n ps and the phase separation density for 
J/t = 1 is np S = 0.745 ± 0.015. 

Emery, Kivelson, and LincJ calculated the phase sepa- 
ration density using the energy per hole, 



where E(Nh) is the total energy of the A^ s -site system 
with Nh holes, and x — Nh/N s is the hole density. In 
Fig. ^| the energy per hole is plotted for all the points cal- 
culated for 42, 50, and 60 electrons. The cubic fits attain 
minima at approximately the same values as the tangent 
constructions. The energies in Fig. ^ are not shifted as 
they are in Fig. ||, and the shell effects are obvious. For 
each number of electrons, the energy is a smooth func- 
tion of the density. However, taking all electron numbers 
together, the energy is a very jagged function: The shell 
effects systematically bias the energies of systems with 
a given number of electrons. Therefore it is essential to 
compare the energies of systems with the same number 
of electrons, thus canceling the unavoidable systematic 
errors. Many previous studies of the 2D t-J model suf- 
fered from shell effects. A different demonstration of shell 
effects is given in Ref. ^l|. 



C. Results in the J < t Region 
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FIG. 10. The ground-state energy per site at J = 0.5t for 
N e = 32 electrons and N a = 36, 49, 56, 64, 72, 81, 90 sites (top 
curve) N e = 42 and iV s = 49,56,64,72,81,90 (second from 
the top), N e = 50 and N s = 56,64,72,81,90 (third from the 
top) and A e = 60 and N a = 49, 56, 64, 72, 81, 90, 100, 110 sites 
(bottom curve), each curve has been shifted downward with 
the respect to the previous by 0.05 In order to distinguish 
them. 

The ground-state energy per site at J = 0.5t for fixed 
number of electrons N e and for various number of sites N s 
can be grouped together. The top curve in Fig. [To| gives 
the energy per site as a function of density n e = N e /N s 
for N e = 32 andN s = 36(6 x 6), 49(7 x 7), 56(7 x 8), 64(8 x 
8) , 72(8 x 9) , 81 (9 x 9) and 90(9 x 10) . The second from the 
top gives the energy per site shifted by a constant amount 
of 0.05 for N e = 42 and A^ s = 49,56,64,72,81,90. The 
third from the top gives the energy per site shifted by an 
0.1 for N e = 50 and N s = 56, 64, 72, 81, 90. The bottom 
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curve gives the energy per site shifted by 0.15 for N e — 60 
and N s = 49, 56, 64, 72, 81, 90, 100(10 x 10), 1 10(11 x 10). 
The unshifted energies are given in Tabic E3. 

TABLE II. The energy per site for J/t — 0.5 and for vari- 



ous 


electrons densities and 


siz6 lclttiC6S. 




N e 


N s 


n e 


E ex trapf -^V s 


Elanc/N s 


32 


64 


0.500 


-1.0102(19) 


-1.0111(18) 




56 


0.571 


-1.0151(24) 


-1.0147(14) 




49 


0.653 


-0.9811(27) 


-0.9802(15) 




42 


0.762 


-0.8890(33) 


-0.8893(32) 




36 


0.889 


-0.7272(11) 


-0.7277(17) 


42 


100 


0.420 


-0.9695(22) 


-0.9698(21) 




81 


0.519 


-1.0158(31) 


-1.0143(22) 




72 


0.583 


-1.0115(50) 


-1.0030(15) 




64 


0.656 


-0.9734(32) 


-0.9734(22) 




56 


0.750 


-0.8918(39) 


-0.8851(5) 




49 


0.857 


-0.7715(24) 


-0.7733(22) 


50 


100 


0.500 


-1.0015(37) 


-0.9925(11) 




90 


0.556 


-0.9975(42) 


-1.0009(29) 




81 


0.617 


-0.9898(54) 


-0.9887(68) 




72 


0.694 


-0.9310(24) 


-0.9340(36) 




64 


0.781 


-0.8594(37) 


-0.8565(59) 




56 


0.893 


-0.7232(22) 


-0.7242(15) 


52 


100 


0.520 


-1.0090(40) 


-1.0102(44) 




90 


0.578 


-1.0020(58) 


-1.0032(80) 




81 


0.642 


-0.9791(46) 


-0.9774(42) 




72 


0.722 


-0.9110(44) 


-0.9107(62) 




64 


0.812 


-0.8158(29) 


-0.8153(20) 




56 


0.929 


-0.6688(31) 


-0.6642(6) 


60 


110 


0.545 


-1.0069(33) 


-1.0041(15) 




100 


0.600 


-0.9929(24) 


-0.9911(20) 




90 


0.667 


-0.9676(32) 


-0.9687(42) 




81 


0.741 


-0.8938(11) 


-0.8940(12) 




72 


0.833 


-0.7947(31) 


-0.7918(18) 




64 


0.938 


-0.6623(25) 


-0.6605(31) 



Here again, we can greatly eliminate the shell effects 
by examining several size lattices but keeping the num- 
ber of electrons fixed. We first fit each curve generated 
for fixed N e with a cubic spline where the Heisenberg 
point has been excluded from the fit. Next we find the 
point (n e = 1, eH{N e )) on the graph, where eH(N e ) is the 
energy per electron for the Heisenberg antiferromagnet 
calculated on a finite-size system with the same number 
of electrons N e (as discussed previously). Next we con- 
struct the tangent to the spline which fits our points. The 
value of the density at which the line is tangent to the 
spline gives us an estimate of the phase separation den- 
sity for this value of J/ 1. These values extracted from the 
different sets of energies which correspond to the same 
number of electrons are given in Table III . Individually, 
each number of electrons is consistent with a value of n ps 
near n e — 0.84 (see Fig. [l0]). Clearly the 42 electron data 
doesn't prove that there is a clear tangent at this value 
of n e , but the data is consistent with this value. There- 
fore we conclude that the finite-size effects are small in 
our method of determining n ps and the phase separation 



TABLE III. The phase separation density at J/t — 0.5 
determined by keeping the electron number fixed and varying 
the lattice size. 

n cr N e 



0.831782 
0.838086 
0.840684 
0.847915 
0.857529 



0.00309639 

0.0136867 

0.00238198 

0.00233589 

0.00913804 



32 
42 
50 
52 
60 



density for J/t = 0.5 is n ps = 0.843 ± 0.015. 
-0.4 
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FIG. 11. The energy per site at J = 
lattice for 32, 42, 50 and 60 electrons. 



0.5t for an 8 x 8 



We wish to demonstrate the significance of shell effects. 
Let us select from our results of Table for J/t = 0.5 
those which correspond to the same size lattice 8x8 for 
N e = 32,42,50,60. They are shown in Fig. [0]. Notice 
that even though these data also give the same phase sep- 
aration density within error bars as that determined by 
our method described before, the shell effects are large. 
Such deviations from a smooth curve could lead to draw- 
ing the wrong conclusions about phase separation bound- 
aries. 



For completeness in Fig. 12 the energy per hole is given 
for all the points calculated for 32, 42, 50 and 60 elec- 
trons. The curve attains a minimum at approximately 
the same value as that determined by the tangent con- 
struction at the cubic polynomial fit of the energy per 
size for a given number of electrons. Notice, again, the 
shell effects. 



D. Results near J/r . 

Boninsegni and ManousakisE! (BM) found a critical 
value J^ ~ 0.27t of J/t below which there is no two-hole 
d-wave bound state. This value of J/ 3 was determined 
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FIG. 12. The energy per hole at J = 0.5t for 32, 42, 50 
and 60 electrons. 



by calculating the binding energy for two holes on lat- 
tices up to 8 x 8. BM noticed that because the bound 
state wave function decays exponentially with distance 
the finite size effects were rather small. They did, how- 
ever, pursued a finite-size analysis from which they de- 
termined Jf . Nevertheless their calculated value of the 
two hole binding energy at J/t = 0.7 was A/t = 0.31(03) 
and at J/t = 0.4 A/t = 0.12(04). Thus, we choose the 
J/t = 0.3 to examine the question of phase separation 
believing that this value is very close to the critical value 
J?. 

In Fig. [ll| we give the ground state energy as a function 
of the electron density for 50,52 and for 60 electrons for 
J/t = 0.3 as three shifted curves. Notice that the values 
of J c /t determined from the these sets of data are very 
close. We obtain: n e = 0.877 ± 0.010. 



TABLE IV. 


The energy per site for J/t — 


0.3 and for var- 


ious electrons densities and 


size lattices. 




N e N a 


n e 


E ex trap/ s 


Elanc/Ng 


50 90 


0.556 


-0.9494(28) 


-0.9469(25) 


81 


0.617 


-0.9199(37) 


-0.9174(23) 


72 


0.694 


-0.8469(25) 


-0.8482(25) 


64 


0.781 


-0.7381(42) 


-0.7419(48) 


56 


0.893 


-0.5510(13) 


-0.5523(15) 


52 90 


0.578 


-0.9407(31) 


-0.9372(38) 


81 


0.642 


-0.9064(41) 


-0.9099(48) 


72 


0.722 


-0.8110(38) 


-0.8095(29) 


64 


0.812 


-0.6881(34) 


-0.6852(49) 


56 


0.929 


-0.4806(18) 


-0.4777(9) 


60 90 


0.667 


-0.8814(11) 


-0.8815(13) 


81 


0.741 


-0.8010(13) 


-0.7976(23) 


72 


0.833 


-0.6538(34) 


-0.6496(13) 


64 


0.938 


-0.4640(23) 


-0.4615(9) 



For J/t = 0.3 in Fig. |lj there is a minimum at n ps = 
0.12 which agrees very well with the value obtained from 
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FIG. 13. The ground-state energy per site at 

J = 0.3t for 50, 52 and 60 electrons and lattices of 
sizes N 3 = 56,64,72,81,90, N s = 56,64,72,81,90 and 
N s = 64, 72, 81, 90 respectively. 



the tangent construction. 

There are no published GFMC results for the two-hole 
case for J/t = 0.3. We can obtain an estimate for the 
single-hole energy by fitting the calculated values for that 
as a function of J/t to a form E — E + aJ 2 / 3 . The 
two-hole binding energy for J/t = 0.3 can be estimated 
using the fownula which were used by Boninsegni and 
ManousakiaLa to obtain the critical value of jf . Thus, 
assuming that holes are bound in pairs and they form 
a dilute gas of hole-pairs, we can obtain a value for the 
energy per hole in such a case. This value of this energy 
is higher than the value of the energy per hole at the 
minimum of our curve in Fig. |TJ. This is another indica- 
tion that there is more binding energy gained due to the 
phase separation of the pairs of holes from the electrons 
in an antiferromagnetically ordered state. 



E. Results below jf . 

Here we examine the situation below the critical value 
jB for two hole d-wave bound state in the 2D t- J model. 
We shall examine the energy at J/t = 0.2. First of all the 
ground state energy per site for N e — 50 and N e = 60 
and for various size lattices is shown in Fig. ^| and is 
given m Table 

Notice again that the values of J c /t determined from 
the two sets of data are very close, we find: n ps = 0.909 ± 
0.008. 

Below jB where there is no two-hole d-wave bound 
state, (assuming there are no bound states in other chan- 
nels) if there is no phase separation the minimum energy 
per hole should be the single-hole energy at zero hole 
density. At J/t = 0.2 the single-hole energy was also cal- 
culated by Boninsegni and Manousakis for an 8 x 8 and 
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TABLE V. The energy per site for J/t = 0.2 and for vari- 
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FIG. 14. The energy per hole at J = 0.3t for 50, 52 and 
60 electrons. 
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FIG. 15. The ground-state energy per site at 

J = 0.2t for 50 and 60 electrons and lattices with sizes 
N s = 56, 64, 72, 81, 90 and N s = 64, 72, 81, 90 respectively. 



10 x 10 size lattices. This value of the energy is shown 
if Fig. |l6| and it is clearly higher than the minimum of 
the energy per hole curve which occurs at approximately 
hole density of x = 0.09. 



V. PHASE DIAGRAM OF THE t-J MODEL 

In Fig. [I?] we present as a function of J/t the minimum 
energy per hole (solid line) (the minimum of the energy 
per hole versus density for a given value of J/t). We also 
plot the single hole energy (obtained by Boninsegni and 
ManousakiaiS) as a function of J/t which is the energy 
per hole in the case of isolated non-interacting holes in 
the system (dashed line). In addition, the energy per hole 
is compared with the esergy per hole obtained by Bonin- 
segni and Manousakist3 from calculation of two holes in 
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electrons densities and 


size lattices. 
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FIG. 16. The energy per hole at J = 0.2t for 50 and 60 
electrons. The single-hole energy as obtained from Boninsegni 
and Manousakis is also plotted for 8x8 and 10 x 10 size 
lattices. 



the t-J model. The energy per hole in this latter calcu- 
lation gives the energy per hole in the case of isolated 
bound hole pairs (dotted line). Notice that while the 
dashed line and the dotted line meet at J/t ~ 0.3, the 
minimum at the phase separation density and the dot- 
ted line do not meet. Notice that the additional energy 
gained to to phase separation decreases with decreasing 
J/t as expected. 

In Fig. |l§| we show the phase separation boundary ob- 
tained for all values of J/t using the method described 
in the present paper and the Maxwell construction. In 
Table VI we give the phase separation boundary as de- 



termined for various values of J/t from the various size 
lattices and number of electrons. 

A more complete phase diagram for the 2D t-J model 
as a function of J/t and doping was given in Fig. 3 of 
Rcf. [lj]. That phase diagram is also accurate in theiow 
density region where exact calculations can be done.EJ 
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FIG. 17. The energy per hole at the density where the 
phase separation minimum occurs as a function of J/t (solid 
line). This is compared to the energy per hole obtained from 
the single hole calculation of Boninsegni and Manousakis ('92) 
(dashed line) and to the energy per hole obtained from the 2 
hole calculation of Boninsegni and Manousakis ('93). 



VI. COMPARISON WITH OTHER 
CALCULATIONS 




FIG. 18. The phase separation boundary as calculated 
using the present method and the Maxwell construction. 



is that the interesting region of J / t is either next to the 
phase separation boundary or inside the phase separated 
region. In both cases phase separation fluctuations could 
play an important role in the mechanism for supercon- 
ductivity in the copper oxides. 



In Fig. our phase diagram is compared to the recent 
fixed node Monte Carlo calculations of Callandra, Becca, 
and Sorella (CBS)tZI and to the high temperature series 
expansion calculations of Putikka, Luchini, and Rice.E-3 
Notice that our phase diagram and that of CBS are very 
close except in the delicate physical region of small J/t. 
Therefore, we can draw a relatively strong conclusion 
from this comparison: The findings drawn from the early 
studies of the t- J model that the physical region of the 
model is safely away from the phase separation boundary 
are not correct. What our work and the work of CBS find 

TABLE VI. The phase separation boundary as calculated 
using the present method. The last value with n = is derived 
analytically in Ref. |l4]. 



J 


n 


a 


0.1 


0.9484 


0.017 


0.2 


0.909 


0.008 


0.3 


0.877 


0.010 


0.5 


0.845 


0.015 


1.0 


0.730 


0.016 


1.25 


0.624 


0.010 


1.5 


0.568 


0.027 


2.0 


0.439 


0.008 


2.5 


0.296 


0.004 


3.0 


0.145 


0.0016 


3.25 


0.0662 


0.0006 


3.4367 
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FIG. 19. Comparison f&f our phase-se«aration boundary 
with that of Putikka et alri and of CBSO 

There is an important difference between our results 
and those of CBS. Our results indicate that phase sepa- 
ration in the t-J model is present for all J/t, while the 
conclusion of CBS is that there is a finite value of J ~ OAt 
below which there is no phase separation. The reason for 
this disagreement is that this region requires a very high 
degree of accuracy in the numerical results. We would 
like to discuss the results of CBS where they claim that 
at J/t — 0.4 there is no phase separation for lattices of 
size N s = 98. In Fig. |fj| we plot the results of CBS for 
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FIG. 20. Quadratic fits of the results obtained by CBsEl 
for systems with 50 sites (solid circles) and 98 site (open 
squares) to a quadratic polynomial. The result for the low- 
est value of x for the 98 site system was not included in the 
original publication by CBS. 

this value of J/t for 50 sites (solid circles) and 98 sites 
(open squares). The result for the lowest value of x for 
the 98 site system was not included in the original pub- 
lication by CBS. CBS were kind enough to calculate it 
at our request and to communicate it to us. Without 
using that point, CBS concluded that the fact that we 
found PS at J/t = 0.4 was a finite-size effect because the 
energy per site in their largest size system had no min- 
imum. With the most recently calculated point for the 
98-site, eh[x) has a minimum at x c ~ 0.072. This is close 
to our value of about x c » 0.1 for J/t = 0.4. 

Let us now examine more specifically the results in 
Fig. ||^. We have labeled by 2, 4, and 6 the points which 
correspond to 2, 4, and 6 holes in the 50 and 98 site 
lattices. Notice that the energy of 4 holes is the same 
within error bars in both lattices. The same is true for 
the 6 hole case. Thus, the energy for 2,4 and 6 holes seems 
to be independent of the size of the lattice within error 
bars. This can be a either a) a genuine characteristic 
of presence of phase separation where the two, four and 
six-hole bubbles in a much larger system do not feel the 
size effects because they are self bound at a characteristic 
size much smaller than the total system or b) a result of 
shell-effects which we have discussed and are minimized 
in our calculation or c) the calculation of CBS has larger 
systematic or statistical errors than those reflected by 
their error bars. 

White and Scalapino (WS) calculated the energy per 
hole on systems with cylindrical boundary conditions, 
that is, systems with open boundaries in. one, direction 
of the lattice and periodic in the otherS'EacZl They es- 
timate the energy per hole, Eq. d53|), by comparing the 
energy of a system with holes to the energy of the same 
system with no holes. In Ref. [l7[ WS argue that their 



approach is more accurate than that obtained by other 
methods simply because it gives a lower energy per hole. 
However, the energy per hole calculated in this way on 
systems with open boundary conditions is not variational 
and, as shown below, can artificially underestimate the 
energy per hole. 

Systems with open boundary conditions can be made 
from fully periodic systems by removing a row of bonds. 
Clearly this process disrupts the periodic ground state 
and raises the energyo Both the energy of the system 
with Nh holes, E(Nf l ), and the energy of the system with 
no holes, £7(0), increase with open boundary conditions, 
but generally not by the same amount. The system with 
holes has more degrees of freedom than the no-hole sys- 
tem, allowing it to respond more effectively to the bro- 
ken bonds. For example, a system with holes has free- 
dom to twist the antiferromagnetic order parameter at 
the boundary required by certain phase separated states. 
An example of such a state is a single "stripe" shown in 
Fig. 551 In this example, a striped or structured phase 
separated state is stabilized in the middle of the system. 

Thus, one expects that the energy per hole obtained 
with open boundary conditions (using as a reference state 
the no hole energy with open boundary conditions) can 
be lower than the exact energy per hole obtained with 
periodic boundary conditions (using as a reference state 
the no hole energy with periodic boundary conditions). 

In Fig. H|, we compare the results of various calcu- 
lations on similar size lattices for the energy per hole 
for 1, 2, 4 holes and at our phase separation minimum. 
The results for one and two holes are taken—from the 
work of Boninsegni and Manousakis (BM)E3'ILj. The fi- 
nite size effects are smaller than the size of the symbols. 
In addition, the tesult for 2 holes for a 50 site cluster 
reported by CBScZl for J/t = 0.4 is shown as an open 
square. Notice the agreement between BM and CBS 
(both used periodic boundary conditions). The value for 
the single hole energy obtained by WS is systematically 
lower than the value for the periodic lattice. The cylin- 
drical boundary conditions used by WS frustrate the no 
hole state, and, as a result, the energy of the no hole 
state obtained by WS is much higher than that used by 
HM and CBS. WS's calculation gives a total energy of 
£7(0) = -35.66 (in units of t and here J/t = 0.5) for the 
no-hole state on the 8x8 lattice, while for a periodic 
8x8 lattice the energy which we (and CBS) use is much 
lower, £7(0) = —37.56. WS's total energy for 4 electrons 
in a 8 x 8 lattice is £7(4) = -41. 028± 0.075, while we find 
£7(4) = -42.23 ± 0.12 on a periodic lattice. Thus WS 
obtain a value for the energy per hole e^(4) = —1.34, 
while our result corresponds to efe(4) = — 1.17 ± 0.03. 

Qualitatively similar conclusions can be drawn for the 
case of a single hole. The results of BM for a single are 
shown by the dashed line in Fig. ^2[ Clearly, WS's single- 
hole energies are below those also. This lowering of the 
energy can only be understood by the frustrating effect 
of open boundary on the antiferromagnetic state. 

Finally notice the very small difference between the 
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FIG. 21. A two-dimensional stripe-type phase separated 
state. State (b) has a pi-phase shift which accommodates the 
hole motion along the stripe but frustrates the AF order in the 
case of periodic boundary conditions. In state (b) this twist 
of the order parameter has no magnetic energy cost with open 
boundary conditions along the x direction. Thus, periodic BC 
conditions in this case frustrate either the hole motion along 
the "stripe" (state (a)) or the antiferromagnetic state (state 
(b)) along the boundary bonds. 
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FIG. 22. Comparison of the energy per hole e h (N h ) at 
Nh = 1, 2, and 4 holes and at the phase separation min- 
imum. The dashed, dotted, and solid lines are polynomial 
fits to e h (N h = 1) from BM, e h (N h = 2) from BM, and e h 
at the phase separation minimum from HM. Notice that be- 
cause of the cylindrical boundary conditions which frustrates 
the no-hole state, WS tend to get more lowering of the en- 
ergy when they introduce holes. For these size lattices the 
finite-size effects on the 1 and 2 hole calculations are smaller 
than the symbol size. For comparison we have also placed the 
result of CBS for 2 holes in a 50 site lattice which is available 
for J/t = 0.4. Notice that the CBS and BM results are nearly 
identical. 



energy per hole in the 2 hole case and in the 4 hole case 
obtained by WS at J/t = 0.35. They find: e h (2) = 
— 1.72 and energy per hole for a stripe —1.737 at the 
optimum doping of 4 holes per stripe. The difference is 
very small and suggests that the WS striped state is only 
a manifestation of frustrated phase separation. 

VII. CONCLUSIONS 

We have developed an efficient Greens Function Monte 
Carlo method for fermions on a lattice that iteratively 
projects out the ground state with no approximations. 
Fermionic minus-sign fluctuations are controlled by using 
all powers of the projection operator up to some maxi- 
mum and extrapolating to infinite power. Starting from 
a good initial state allows us to converge before the sta- 
tistical errors become too large. This technique comes 
also with solutions to a number of other technical prob- 
lems such as a) enabling the guided random walk to walk 
through the nodes with an 0(N 2 ) algorithm using the 
idea of a "detour walk" b) using a single walker to com- 
pute all the desired powers of the projection operator 
(H — W) n \ where m = 0, 1, ...,p mal , simultaneously. 

This technique is applied to the two-dimensional t-J 
model to investigate its phase diagram. It is found, con- 



trary to many previous studies that there is phase sepa- 
ration (PS) at all interaction strengths of the t-J model. 
The signal for phase separation is clear when one over- 
comes the following difficulties: 

First, the Maxwell construction is the cleanest and 
strongest signal for PS because it suffers the least from 
finite-size effects. Second the shell-effects can mask the 
signal because the energy as a function of density of dif- 
ferent numbers of electrons on a fixed lattice is not a 
smooth curve. The kinetic energy jumps discontinuously 
as electrons are added to successive shells. Therefore, 
we have chosen to keep the electron number fixed at a 
closed shell configuration and to change the size of the 
lattice. The number of electrons which form closed shell 
configurations depends on the boundary conditions. To 
generate as many as possible "magic numbers" of closed 
shell configurations we have used four types of bound- 
ary conditions. Periodic with or tt phase shifts at the 
boundary in each of the x and y directions. 

We find that for any value of J/t the energy per site 
e(n e ) as a function of electron density n e for finite size lat- 
tices does not remain a concave function at high electron 
density. There is a value of the density n ps (J/t) where 
a straight line starting from the no-hole energy per site 
is tangent to the curve e(n e ) at n e = n ps . While the 
energy e(n e < n ps ) does not change significantly with 
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system size, the energy of a finite system in the phase 
separated regime e(n e > n ps ) changes with system size 
and approaches this tangent line in the infinite size limit. 
We interpret this as evidence for phase separation at all 
values of J/t. The fact that the function e(n e ) does not 
remain concave in our calculation above n ps (J/t) can be 
explained by the energy cost of forming an interface be- 
tween the two phases in our finite system. 

Our results have been compared to the most works of 
Calandra, Becca, and Sorellaj23 and we find very close 
agreement. These comparisons indicate that the early 
conclusions that the critical J c /t for phase separation is 
far away from the physical value of J/t are largely invalid. 
This comparison also indicates that J c /t is very small 
and may vanish. Wej-discuss recent comparison by White 
and Scalapino (WS)cZl of our numerical results to theirs. 
In that comparison, WS use the variational principle to 
argue that their results are more accurate because the 
energy per hole in lower. However, we demonstrate that 
one should expect the exact energy per hole on periodic 
lattices to be higher than that obtained with the cylin- 
drical boundary conditions used by WS. Thus on such 
different systems, a lower energy cannot be used as a cri- 
terion for the accuracy of an approach. In addition, we 
interpret the results of WS as evidence for phase separa- 
tion and the appearance of stripes in the t-J model as a 
finite-size effect. 
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APPENDIX A: INVERSE UPDATE THROUGH 
NODES 

To evaluate the determinant in the trial state, we 
use the usual "inverse update" trick first applied to 
condensed matter systems by Ceperley, Chester, and 
KalosJiJ We calculate the determinant and inverse of the 
matrix ( p6|) together at the start of each run, an opera- 
tion taking 0(N 3 ) steps for aJVxJV determinant. Then 
with each single particle move in the random walk, we 
update the determinant in O(N) steps and the inverse in 
0(N 2 ) steps. 

Starting with the matrix D and its inverse I, suppose 
we change row I of the matrix to Dij — > rj. Since the 
inverse is the transpose of the matrix of cofactors nor- 
malized by the determinant, 



h 3 =cof ji (D)/|D|, (Al) 

the ratio of the determinant of D before and after the 
change is 

* s w = 5> J *- (A2) 

3 

The new inverse matrix is given by 

I'ij = h (l + - ]T r k I kj , (A3) 

and one can easily confirm . ^'ij^'jk = ^ifc- Changing 
one column of the matrix results in a similar update for 
the inverse. 

The algorithm is straight forward, and has been used 
in many GFMC studies in the continuum and in Varia- 
tional Monte Carlo on a lattice. However, it cannot be 
used directly with GFMC on a lattice since the random 
walk steps directly on nodes for a significant fraction of 
steps. When the matrix becomes singular, its inverse is 
undefined, and the algorithm breaks down. 

One way around this problem is to recalculate the 
determinant and inverse after walking through a node. 
However, in a reasonably dense system, a large faction of 
steps will land on nodes, and the running time will scale 
as 0(N 3 ). 

We developed a new 0(N 2 ) technique to hop 
over nodes without recalculation of the determi- 
nant or inverse. The essence of the method is 
this: Let us suppose that the random walk vis- 
its a node; namely, the particles were in a config- 
uration R = (ri T ,r 2T , .:,r N /2i,r u ,r 2 i, fjv/24.) and 
by moving a particle, say the first up-spin particle 
from position r\-\ to the determinant defined by 

Eq. (|2^) is zero for the new configuration R' = 
{7 lv r2^...,r N/ 2i,r ll ,r 2 i,...,r N/ 2i). That is a problem 
for the application of the inverse update. In order to 
move to the next configuration, say where particle 2 is 
positioned at f 2 ^ and this corresponds to a new configu- 
ration R" = (f lv f 2V ...,r N/ 2 1 ,r ll ,r2i,...,fNi), we need 
the inverse matrix I = I(R') for the configuration R' and 
this does not exist because the determinant D(R') = 0. 
The non-existence of the inverse is no problem for the 
physics because all we need for computing the observables 
is the determinant, not the inverse; the inverse matrix is 
only a tool which saves us from having to recalculate the 
full determinant at each step. However, in a reasonably 
dense system, a large faction of steps will land on nodes, 
and the running time will scale as 0(N 3 ). We have been 
able to use the inverse update technique by making a 
"detour" around the node as follows. The real motion 
of the random walk was R — > R' — > R" and because 
D(R') = we cannot update the inverse to find I(R') 
which we need in order to calculate D(R"). However, 
all we need is D(R") and I(R") independently of how 
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the random walk got there. Let us consider the configu- 
ration R' 2 = (r lh r , 2V ...,r N /2T,r ll ,r 2l ,...,rNi), which is 
obtained from R by imagining that we moved (without 
actually doing it) particle 2 with spin-up to 7^ and let us 

assume that D(R' 2 ) 7^ 0. Since R" can be obtained from 
R! 2 by moving particle 1 to r[p D(R") and I(R") can be 
obtained by imagining that the walk went through R 2 to 
get to R". Thus, in order to calculate D(R") and I(R") 
we only need to calculate D(R' 2 ) and I(R' 2 ) which both 
exist. 

When the random walk generated by the guiding func- 
tion hits a state or series of states where the determinant 
of the trial function vanishes, we generate a "detour" 
walk around the region where the matrix is singular, re- 
joining the guiding walk when the determinant is non- 
zero again. 

To choose the detour walk, we simply delay any move 
causing the determinant to vanish and place the particle 
number and its future site at the beginning of a list of 
moves to make. For any subsequent move of a particle of 
the same spin, we try to move the first particle in the list 
to that site. If that move yields a non-zero determinant, 
we accept it and attempt to move the next particle in 
the list in the same manner. We repeat the process until 
either all moves give zero determinant or the list is empty, 
in which case the true determinant is not zero. 

Obviously, the procedure will not produce a non-zero 
determinant when the true determinant is zero. How- 
ever, it is important to prove that the detour rejoins the 
guiding walk at the first step with non-zero determinant. 

We represent the rows of the matrix ( pif ) by D = 
{|ri), |r 2 ), . . . , |r„)} where \fi) represents the row 



For a Fermi liquid state, fl2q ) may be expanded into 
the product of two Slater determinants, and this algo- 
rithm suffices as it stands. However with a pairing trial 
state, this decomposition is not possible, we must con- 
sider moves of opposite spin electrons causing the deter- 
minant to vanish. 

Suppose we find moving either the first up particle, 
changing the first row to D±j — > rj, or the first down 
particle, changing the first column to Dn — > Cj, results 
in a zero determinant for matrix (p6|). 

If In / we move the first row with the first element 
shifted by 1/In, noting 



ID'I = 



n 



i 

in 

D 2 i D 22 



r-2 



Dnl D n2 



IDI 



(A7) 



We then try to change the first column in the standard 
manner. We have artificially changed the upper left el- 
ement of the determinant, but since this element will 
be changed again before we finish the detour walk, the 
change will not affect the true determinant. 

If In = 0, this modified step is no longer possible, and 
we need to prove that the true determinant, 



ID"! 



x r 2 

C2 D 2 2 



D 2n 



D 



u2 



(A8) 



vanishes. Here x is the upper left element after both 
moves. 

We know there exist coefficients a 2 , 0.3, . . . , a n such 



\n) = (a(n-f - rti),a(ft\ - raj.), a{n] - r N/2i )) (A4) that rj = J2 t >2 a i d v for au 3- Since the inverse is related 



which is labeled by fi. Suppose moving the first up par- 
ticle to a new site, changing the first row label to rj — > s, 
yields a zero determinant. Then 



to the matrix of cofactors by (Al), cofn(D) = 0, and 
there are other coefficients 0%, 03 , ■ ■ ■ , Pn such that = 
•£ z > 2 fadij for all j > 2. If J2 t > 2 ^ = ^ th en |D"| = 
trivially. Otherwise let 7; = a, + A/% where 



|s) = a 2 \r 2 ) +a 3 |r 3 ) 



(A5) 



A = 



for some coefficients a2, 0:3, . . . , a n . Let the next random 
walk step move the second particle, changing row r 2 — > t. 
Simply by checking if the matrix D' = {|t), |ra), . . . , |r„)} 
has zero determinant, we can determine if the true matrix 
D" = {|s), |t), Irg), . . . , |r„)} is singular. If |D'| ^ 0, we 
accept the move, swap the particles, and try to move the 
first particle again. If |D'| = 0, 



E 



OLiCi 



(A9) 



and 



T,i>2TA] for a11 3 > 2, 



|t) = &|r 2 )+/%|r 3 ) + ---+/? n |r n ) 



(A6) 



Then x = J2i_ 
so |D"| = 0. 

Again, this argument can be extended to any number 
of delayed moves. By combining the two types of moves 
described in this section, we are able to keep track of the 
true determinant without recalculating the inverse from 
scratch. 



for certain coefficients @ 2 , /?3, ■ ■ ■ , n - Combining (A5) 
and (A6) to eliminate r 2 , we see that |D'| = implies 
|D"| = 0. Thus by simply checking single particle moves, 
we can verify that the determinant of the matrix two 
steps away is zero. The argument is easily generalized to 
any number of delayed moves. 



APPENDIX B: O(N) CALCULATION OF 
SUPEREXCHANGE 

For a determinantal function, the kinetic terms in 
require O(N) steps per particle, so it scales as 0(N 2 



for 



20 



the system. The superexchange term in the t-J model, 
^i~S~ , exchanges two particles, changing both a 
row and a column of the determinant (^6|). In this sec- 
tion, we show how the amplitude of swapping two parti- 
cles may be calculated in O(N) steps. 

Suppose we swap the m'th up electron with the n'th 
down electron. We will modify both row m and column 
n in the determinant. We write the new elements as 
D m j — > Tj and Di n — > Cj. Naturally, r n — c m . One can 
show the ratio of the determinant before and after the 
swap is 



ID 



(Bl) 



rilijCj takes 



Direct evaluation of the sum S = ^ 
0(N 2 ) per pair of neighboring particles. For this reason, 
many researchers evaluate the superexchange term only 
every N Monte Carlo steps E3 

Our trick is to evaluate S once when a pair of particles 
become nearest neighbors, and then to update it in O(N) 
steps for any move not disrupting the pair. 

Suppose the Z'th up electron moves (I ^ m), altering 



so D tj - 
t) and ci 



Sj. The inverse 



c'i takes a new 



row I in the determinant (|2(: 
/ is updated according to 
value. 

We can write the new sum S' in terms of the old sum 
and extra factors as 



n (kiil + -5ij) - -la E s ^ ) c 'j ( B2 ) 



S- 



where jj = J^k s klkj is used in the inverse update. This 
calculation requires only 0(N) steps, so the local su- 
perexchange energy of the system may be evaluated in 
0{N 2 ) time. 



APPENDIX C: ENERGY AT THE PHASE 
SEPARATION BOUNDARY 

In the phase separated state, the t-J model separates 
into two phases, one with all electrons (no holes) and the 
other with some electrons and some holes. The transition 
is continuous: As J is increased in the phase separated 
regime, the electron density in the low-electron-density 
phase decreases while the proportion of Heisenberg phase 
increases. The energy in the partially phase separated 
regime is simply the weighted sum of the two constituent 



energies. Specifically, the energy of the phase separated 
state is given by 



1 — n 



E ps (n, J) = -—^E u (n ps , J) + " ps E H J (CI) 

-L ^r).9 -L Tips 



where E u (n, J) is the energy of the uniform density phase 
as a function of electron density and interaction strength, 
EhJ is the energy of the Heisenberg phase, and n ps (J) 
is the density of the onset of phase separation. 

Across the phase separation boundary, the energy is 
continuous as is its first derivative with respect to density. 
Using this fact, we can show that the derivative of the 
energy in the phase separated regime with respect to J 
is given by 



dE ps (n,J) _ 1-n dE u (n ps ,J) n - n p6 



dJ 



l-n % 



dJ 



l-n. 



E H (C2) 



so the first derivative of the energy with respect to J is 
continuous at the phase separation point, n = n ps (J). 
All terms of the form d»p»0-0 are cance i ec | from this ex- 
pression. Note that for J > J c , where J c is the crit- 
ical interaction strength for complete phase separation, 
n ps (J > J c ) = and E u (n = 0, J) = 0. 
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